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We study numerically the fully nonlinear spherically-symmetric collapse of a self- 
gravitating, minimally-coupled, masslcss scalar field. Our numerical code is based on 
double-null coordinates and on free evolution of the metric functions and the scalar 
field. The numerical code is stable and second-order accurate. We use this code to 
study the late-time asymptotic behavior at fixed r (outside the black hole), along the 
event horizon, and along future null infinity. In all three asymptotic regions we find 
that, after the decay of the quasi-normal modes, the perturbations are dominated by 
inverse power-law tails. The corresponding power indices agree with the integer values 
predicted by linearized theory. We also study the case of a charged black hole nonlin- 
early perturbed by a (neutral) self-gravitating scalar field, and find the same type of 
behavior — i.e., quasi- normal modes followed by inverse power-law tails, with the same 
indices as in the uncharged case. 



1 Introduction 

Until recently, the late-time evolution of non-spherical gravitational collapse was 
investigated primarily in the context of linear theory. The late-time behavior of 
such perturbations has been studied for three different asymptotic regions: (a) at 
fixed r (outside the black hole), (b) along null infinity, and (c) along the future 
event horizon (EH). Qualitatively, the evolution of the linearized perturbations is 
similar in these three asymptotic regions: During the first stage, the perturbations' 
shape depends strongly on the shape of the initial data. This stage is followed by 
the stage of quasi-normal (QN) ringing, Finally, there are also 'tails', characterized 
by an inverse power-law decay. 

gase (a) was first studied by Price El, and cases (b) and (c) by Gundlach et 
It was found that after the QN ringings die out, the perturbations in (a) 
decay according to f-( 2 '+A'+i) ) where /i = 1 if there were an initial static mode, 
and fj, = 2 otherwise. Here, I is the multipole moment of the mode in question, 
and t is the Schwarzschild time coordinate. Tails in (b) decay according to u e 
and in (c) according to Ve^ 2l+fl+1 \ where u e and v e are the outgoing and ingoing 
Eddington-Finkelstein coordinates, correspondingly. 

The numerical simulation of the fully- nonlinear gravitational collapse of a spher- 
ical self- gravitating scalar field was recently, carried out, and the QN ringing and 
the 'tails' were demonstarated for cases (a) Em and (c) □. 

In what follows we shall briefly describe a recent numerical simulation of the 
fully-nonlinear spherical collapse of a minimally-coupled massless scalar field. We 
study all three cases (a), (b) and (c), and obtain values for the power-law indices 
significantly closer to the linear analysis predictions than all previous nonlinear 
simulations. Section describes our numerical method, and Section describes our 
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main results. Further details are given elsewhere 



2 The numerical method 

Our numerical code is based on free evolution of the dynamical equations in double- 
null (DN) coordinates. The constraint equations are imposed only on the initial 
characteristic hypersurface, and just monitored during the evolution. The main 
advantages and properties of the DN approach are: The DN coordinates are very 
well adopted to the hyperbolic character of the field equations, the interpretation of 
the causal structure of the spacetime is trivial, and DN coordinates can be chosen 
such that the metric is regular at the EH. Evolution is obtained by straightforward 
marching along both coordinates. Because in the dynamical equations all second 
derivatives are mixed, a second-order code can be obtained from just two adja- 
cent grid-points in each direction. Therefore, only three grid points are needed in 
each computational cell, and the second-order accuracy is obtained by standard 
'predictor-corrector' technique. The inevitability of some sort of Dynamical Mesh 
Refinement for our code can be demonstrated for Schwarzschild (for any spacetime 
with an EH the reasonings are similar). First, the EH must be included in the 
integration. If it were not, then very quickly the integration domain will get very 
far from the EH, and none of the cases (a), (b) or (c) would be susceptible of ac- 
curate analysis. However, whenever the EH is included in the integration we face 
a fundamental difficulty: As one moved along an outgoing null ray near the EH, 
r tU grows exponentially, where r is the area coordinate. Hence, a small error in the 
location of the EH would quickly blow off. The solution is to make the grid dense 
near the EH (to allow for an accurate integration), but dilute far from the EH (to 
save computation time). This is done by the following algorithms. Point Splitting: 
We check the variation of the metric functions and the scalar field along ingoing rays 
between two adjecent grid points. If the variation is greater than some threshold 
value, we introduce an intermediate grid-point, at which we evaluate the fields by 
interpolation. Chopping: The domain of integration includes the EH. Consequently, 
(with vanishing charge) any fixed Wfmai would lead to a crash into the singularity 
after a finite lapse of time. We solve the problem by simply chopping ingoing rays 
immediately after the apparent horizon, which we find locally. Chopping enables 
the code to run forever while including the EH in the integration. Due to the Point 
Splitting procedure, we obtain a denser grid than necessary far from the EH. We 
introduce Point Removal to remove unnecessary grid-points, by a method which 
essentially is the inverse of Point Splitting. However, in any Kruskal-like gauge g uv 
grows rapidly along ingoing rays, anf therefore Point Removal would be ineffective. 
We cure the problem by introducing Gauge Correction, which changes the gauge to 
a gauge in which g uv changes only slowly. The combination of Point Removal and 
Gauge Correction enables a very effective saving of computation time. Our code was 
checked by the following methods: We compared the results obtained with different 
grid-parameters. We monitored the discrepancy in the two constraint equations. 
We compared the mass function obtained from local differentiation and from evo- 
lution of the 'wave equation' for the mass function. We numerically reproduced 
known exact solutions [Schwarzschild, Reissner-Nordstrom (RN), the homothetic 
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Table 1: Values for the Local Power Indices for cases (a), (b) and (c) in configurations (1) (un- 
charged case) and (2) (charged case). The value in the brackets is the linearized theory prediction. 





(a) : along the EH 


(b): along null infinity 


(c): along r — const 


Config. (1) 


2.98 ±0.01 (3) 


2.002 ±0.003 (2) 


2.99 ±0.02 (3) 


Config. (2) 


2.99 ±0.01 (3) 


1.996 ±0.001 (2) 


2.99 ±0.02 (3) 



Roberts solution □] . All checks indicated stability, convergence and second-order 
accuracy. 

3 Results for the decay of 'tails' 

We used our numerical code for the following two configurations. (1): Collapse 
over Minkowski, which leads to the formation of a Schwarzschild-like black hole. 
(2): Collapse over a unity mass pre-existing charged background (a RN geometry) 
with various values of the charge. We took an initially ingoing pulse of squared-sine 
shape of compact support, and chose the amplitude to be high, such that the final 
black hole mass was Mf « 3.5 in (1) and Mf « 3.8 in (2). We then evolved the 
fields, and probed their values for cases (a), (b), and (c) in both configurations (1) 
and (2). In order to approximate null infinity [for case (b)] we probed the fields on 
an ingoing null ray at v e = 10 6 Mf. We introduce the notion of the local power 
index (LPI). Namely, we calculate the evolution in p = —v (ln<I>)^ for case (c) [v 
should be replaced by t and u for cases (a) and (b), respectively], where $ is the 
scalar field. The motivation behind the LPI is as follows: The power-law behavior 
is just the leading-order term in an expansion in v^ 1 . For any finite value of v there 
will also be nonvanishing contributions of the higher-order terms, which will cause 
a deviation from the integer power index. In addition, any averaging or best-fit 
technique would yield a result which depends on the interval. The LPI would thus 
be a fractional number, but would asymptotically approach the integer value of the 
leading-order power index. Table [l] summarizes our results. The agreement with 
the linearized-theory predictions is remarkable, and better than the results of all 
previous nonlinear analyses. 
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